Accretion and destruction of planetesimals in turbulent disks 
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q ! ABSTRACT 

We study the conditions for collisions between planetesimals to be accretional 
or disruptive in turbulent disks, through analytical arguments based on fluid 
dynamical simulations and orbital integrations. In turbulent disks, the velocity 
dispersion of planetesimals is pumped up by random gravitational perturbations 
from density fluctuations of the disk gas. When the velocity dispersion is larger 
than the planetesimals' surface escape velocity, collisions between planetesimals 
do not result in accretion, and may even lead to their destruction. In disks with 
a surface density equal to that of the "minimum mass solar nebula" and with 
nominal MRI turbulence, we find that accretion proceeds only for planetesimals 
with sizes above ~ 300 km at 1AU and ~ 1000 km at 5AU. We find that accretion 
is facilitated in disks with smaller masses. However, at 5AU and for nominal 
turbulence strength, km-sized planetesimals are in a highly erosive regime even 
for a disk mass as small as a fraction of the mass of Jupiter. The existence 
of giant planets implies that either turbulence was weaker than calculated by 
standard MRI models or some mechanism was capable of producing Ceres-mass 
planetesimals in very short timescales. In any case, our results show that in the 
presence of turbulence planetesimal accretion is most difficult in massive disks 
and at large orbital distances. 

Subject headings: solar system: formation -- planets and satellites: formation 
— accretion, accretion disks — turbulence 
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Introduction 



It is often considered that the evolution of protoplanetary disks and the consequent 
accretion of gas by the central pro tostar are driven by tur b ulent viscosity due to a Ma gneto- 



Rotat ional-Instability (MRI) (e.g. jBalbus fc Hawleylll99ll ). lLaughlin et al.l (120041 ) and lNelson fc Papaloizp- 



( 120041 ) carried out fluid dynamical simulations of MRI and found that the random torques 
due to the tur bulent density flu c tuatio ns give rise to a random walk in semimajor axes of 
planetesimals. iRice fc Armitagd (120031 ) pointed out through model calculations that the ran- 
dom walk expands the effective feeding zone of protoplanets, and may lead to rapid form ation 
of large cores for gas giants. Through a Fokker-Planck treatment, iJohnson et al.l (120061 ) also 
pointed out the importance of the random walk in planet accretion . Ado p ting the semi 



analyt ical formula for the random torque derived by lLaughlin et al.l (120041 ) , lOgihara et al. 



( 120071 ) performed N-body simulations for the late stages of terrestrial planet accretion with 
a disk significantly depleted in gas, starting from Mars-mass protoplanets. They found that 
the MRI turbulence indeed helps to reduce the number of accreted terrestrial planets which 
is otherwise too large compared to our Solar System. 



However, iNelsonl (120051 ) found through direct integrations of t he orbits of p r otopl anets 
in a MRI turbulent disk that orbital eccentricities are also excited. iBritsch et al.l (120081 ) also 
found a similar feature in a self-gravitating disk. While the random walk itself is favorable 
to the growth of protoplane ts by avoiding iso l ation , the excitation of their e ccentricities, 
which had been neglected in IRice fc Armitagd (120031 ) and IJohnson et al.l (120061 ). is a threat 
for planetesim al accretion processes because of increased collision velocity. Unfortunately, 
Nelson! (J2005|)'s orbital integrations were limited to 100-150 Keplerian times and neglected 
collision processes, so that it is not possible to conclude from that work whether planetesimals 
should grow or be eroded in the presence of turbulence. 

In the present article, we explore by which paths planetesimals may have grown to 
planet-sized bodies in turbulent disks. Because the level of density fluctuations due to the 
MRI turbulence is not well determined, we choose to study the qualitative effects of the 
turbulence on the accretion of planetesimals and their dependence on the key parameters 
of the problem, in particular the progressive removal of the gas disk. In §2, we summarize 
the conditions for the accretion and destruction of planetesimals in terms of their orbital 
eccentricities. In §3, we analytically derive the equilibrium eccentricity for which the exci- 
tation due to turbulence is balanced by damping due to tidal interactions with the disk gas, 
aerodynamic gas drag, and collisions. Comparing the equilibrium eccentricities with critical 
eccentricities for accretion and destruction, we derive critical physical radii and masses of 
planetesimals for accretion or destruction. The results are applied to viscously evolving disks 
(§4). We then discuss possible solutions for the problem of the formation of planetesimals 



and planets (§5). 
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Accretion and destruction conditions 



We summarize the accretion and destruction conditions below. From energy conserva- 
tion, the collision velocity (i> co ii) between two planetesimals (labeled 1 and 2) satisfies 
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where Mj and Rj are the mass and physical radius of a planetesimal j (j = 1, 2), and t> re i is 
their relative velocity when they are apart from each other. When the velocity dispersion of 
planetesimals v^ ap is larger than the Hill velocity that is given by (M/3M») 1/,3 uk, where M* 
is th e mass of the host star and Vk is the Keplerian velocity, v re \ is approximated by t> disp 

(e.g., T ' : """'" " """""" 



Ida fc Nakazawalll989l ; Ildalll990l ). The total energy then becomes: 
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where v esc is the (two-body) surface escape velocity defined by 
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Collisional dissipation decreases the energy by some fraction of v^ oll /2. If f disp <ti v esc , 
the collisional dissipation results in E < after collision. On the other hand, E is likely 
to be still positive after a collision with t>di sp ^> f esc . T hus, for moder ate dissipation, the 
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condition for an accretional collision is t>disp 

eccentricity e ~ v disp /vK, a collision should result in accretion for e < e 



Ohtsuki Ill993l ). Since the orbital 
where 
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In the above relation, p p is the bulk density of the planetesimals and M ~ Mj (for simplicity, 
Mi ~ M 2 is assumed). The physical radius R is given by 

7.8 x 10 8 (M/M e ) 1/3 (p p /3gcnr 3 ^ 1/3 
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A collision results in destruction if the collision velocity is such that the specific kinetic 
energy of a collision (v 2 ol /2) exceeds 
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wher e Qo is the material strength, B ~ 0.3-2.1, a p± —0.4, and 6 — 1-3 (Benz fc Asphaug 
19991 ). For basalt rocks or water ice, Qo — 10 7 -10 8 ( IBenz fc Asphauglll999l ). but it can take 
a significantly smaller value for loose aggregates, lj 

We adopt Qq = 3 x 10 7 as a nominal value. Self-gravity (the second term in the r.h.s.) 
dominates the material strength when R > 100 m. In this regime, adopting B ~ 1, a collision 
results in destruction for e > edis, where 
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3. Equilibrium eccentricities 



We first derive the equilibrium eccentricities of planetesimals at which the excitation by 
the MRI turbulence is balanced by damping due to drag and/or collisions. Comparing the 
estimated eccentricities with e acc and e^s, we then evaluate the outcome of collisions between 
planetesimals as a function of planetesimal size, turbulent strength, and surface density of 
disk gas. 

For an easy interpretation, we provide in this section analytical relations based on a 
model in which the gas and solid components of disk surface density are scaled with the 
multiplicative factors f g and fd- 
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where 77 ice ~ 3-4 is an enhancement factor of S^ due to ice condensati on. If f g = fd = 1, S 9 
and Erf are 1.4 times those of the minimum mass solar nebula model (JHayashilll98ll ). 



In thi s section, we a lso use the disk temperature distribution obtained in the optically 
thin limit (JHayashilll98lf ). 
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-'^For porous materials, Qo is rather higher (W. Benz, private communication). 
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where L* and L & are the stellar and solar luminosities, respectively. The corresponding 
sound velocity is 



1.1 x 10 5 



cm/s. 



(11) 



Since the disk scale height is given by h = v2c s /$7k (assuming that T is vertically uniform 
in the disk), eqs. (JED and fill I) yield the disk gas density at the midplane as 



-11/4 



Pg 



V^h 



2 x 10- 9 /, (— j gem 



(12) 



3.1. Excitation 



The orbital eccentricities of planetesimals are pumped up both by the random gravita- 
tional perturbations from density fluctuations of disk gas, as well as by mutual gravitational 
scattering among planetesimals. Assuming planetesimals have equal masses, the ir orbital 
eccen tricities should be excited to at most ~ e acc by the mutual scattering (e.g., ISafronov 
19691 ) . As will be shown below, the value of this excentricity is smaller than that due to 
the turbulent excitation, except for very large planetesimals (~ 10 3 km or larger), and/or in 
the case of significantly depleted gas disks. For simplicity, in this work we choose to neglect 
the possibility that mutual scattering dominates over turbulent excitation. It is therefore 
important to note that our results may be slightly optimistic when concerning the possibility 
of accretion of massive planetesimals. 

The orbital eccentricities that resu lt from the turbulent density fluctuations in the disk 



are provid ed by lOgihara et al. 



formula by lLaughlin et all (120041 ) , as 



(120071 ) on the basis of orbital integrations with empirical 
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where E 9) i is S 9 at 1AU with f g = 1 (e g. [8]) and 7 i s a no n-dimensional parameter to 



1 — 1 r — ^ ■ — 1 1 1 ■*■ 

express the disk turbulence. Although lOgihara et al.l (120071 ) showed the results only at 



2 Although lOgihara et al.1 (J2007I ) suggested that eq. (fl3l) may be enhanced by a factor 10 by the inclusion 
of to = 1 modes, the to = 1 modes actually enhance only the amplitude of random walk in semimajor axis 
(Aa) but not the eccentricity. Since higher to modes fluctuate over shorter timescales, they tend to cancel 
out on the orbital period of a planetesimal. For these modes, Aa/a, which is due to time variation of the 
potential, is much smaller than Ae, because the latter is also excited by the non-axisymmetric structure. 
The inclusion of slowly varying to = 1 modes enhanc e s Aa/ a up to the order of ~ e. On the other hand, the 
definition of T in eqs. (5) and (34) of lOgihara et al.l (|2007l ) should be multiplied by tt. We use eq. (Til?]) for 
the eccentricity excitation, which is consistent with an orbital calculation including to = 1 modes (figure [1} . 
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~ 1AU, we here added a dependence on r using scaling arguments (see the Appendix). 
Orbital integration for other r show a consistent dependence. Note that Aa/a ~ e, where 
Aa is the amplitude of random walk in semimajor axis. Since e<l, the radial distance r 
and semimajor axis a are identified here. 



From the simulation results by lLaughlin et al.l (12004 ). the value of 7 may be ~ (l/3)(Sp/p) 
10~ 3 -10 -2 for MRI turbulence. In this paper, we use 7 = 10~ 3 as a fiducial value. Inter- 



estingly, with a quite different approach, iJohnson et al.l (120061 ) derived a similar formula for 
Aa/a with the same dependences on r, E 9 and t. If e ~ Aa/a, their formula is consistent 
with ours. They suggested that 7 ~ a or a l / 2 h/a where h i s disk scale height and a is the 
parameter for the alpha prescription for turbulent viscosity (jShakura &: Sunyaevlll973l ). For 
a = 10~ 3 -10~ 2 , their estimate is also similar to our fiducial value. 

The top panels in fig. [T] show the results of an orbital integration with 4-th order Her- 
mite scheme for the evolution of e and Aa/a with turbulent perturbations but without any 
damping. Five independent ru ns with different ran dom number seeds for the generation of 
turbulent density fluctuations (jOgihara et al.l 120071 ) are plotted in each panel. The initial e 



and i are 10~ 6 . For f g — 1, 7 = 0.01, and r = 1AU, as used in fig. [TJ eq. (TT31) is reduced to 
e ~ 10~ 3 (t/lyear) 1 / 2 . To highlight the effect of turbulence, we used a larger value of 7 than 
the fiducial value. The evolution of the root mean squares of the five runs in fig. [1] agrees 
with eq. (fi~3l) within a factor of ~ 2. From eq. (fi~3l) . the excitation timescale is 
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years. 
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3.2. Damping 

The eccentricity damping processes are i) tidal interaction with disk gas, ii) aerody- 
namical gas drag, and iii) inelastic collisions. The tidal damping timescale (i) is derived by 
Tanaka fc WardPhooi ) as 
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The gas drag damping timescale (ii) is derived by lAdachi et al.l (119761 ) as 
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For simplicity, we evaluate the damping timescale due to inelastic collision as the mean 
collision time of planetesimals, assuming that all the planetesimals have the same mass M. 
Since in the size distribution caused by collision cascade, collisions with comparable-sized 
bodies and those with smaller ones contribute similarly, the neglection of the size distribution 
may not be too problematic. Since we look for the conditions in which collisions are non- 
accretional, we consider the case with v d i SV > v csc . Assuming that the gravitational focusing 
factor [1 + (v csc /vdis P ) 2 } ~ 1, the collision damping timescale is 

Wi-^ ^ 2 x lO 7 /^ 1 ( ^f) ( T ^) (^-Vyears, (17) 

nirR 2 v disp \M e J V 3 S cm / vlAU/ 

where n is the spatial number density of planetesimals. Note that n ~ (T>d/ M) / (v disp /Q K ) . 

3.3. Equilibrium eccentricity 

We now equate eq. (fH|) with eqs. ffT5|) . (fT6j) . and (IT7|) . respectively, to obtain an equilib- 
rium eccentricity for each damping process. For simplicity and to a good approximation, the 
actual equilibrium eccentricity can be approximated as the minimum of the three equilibrium 
eccentricities. From eqs. (fill) and (fl~5l) . 
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With eq. (PT6D, 
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For / g = 1, 7 = 0.01, and r = 1AU, eq. (JTHJ) predicts that e drag ~ 0.045 for M/M e = 10~ 6 
and e^rag — 0-01 f° r M/M 9 = 10~ 12 . An orbital integration in the middle and bottom panels 
in fig. [I] shows that the results agree with the analytical estimate within a factor ~ 1.5. With 
eq. ((ED, 

e coll * 3.2 x 10V,(/rfk.)-^ 7 (l^) 1/6 (siT^) 173 (iiu) 5/4 

* 3.6/ p (/^)-^ 7 (^) 1/2 (si^) 576 (tm) 574 - 

In figure [2], the equilibrium eccentricity, e eq = min(e t idai, edrag, e co ii) , is plotted with 
solid lines as a function of the planetesimal radius R, the corresponding planetesimal mass 



being M — 2.1 x 10~ 3 (i?/10 3 km) 3 (p p /3gcm~ 3 )M e , Note again that the effect of mutual 
planetesimal scattering is neglected. For bodies with more than Lunar to Mars masses, tidal 
damping is dominant. This yields a decrease in the equilibrium eccentricity with increasing 
planetesimal radius for R > 100km. For smaller mass bodies, gas drag damping dominates 
tidal damping and the equilibrium eccentricity increases with increasing R. For the smallest 
planetesimal sizes (the regions with the slightly steeper positive gradient), collision damping 
is dominant, but with a significant contribution of gas drag damping. 



The limiting mass and radius at which ed rag (eq. pj5]) and e acc (eq. [I]) cross are 

-V4 , . 9/8 
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The accretion of planetesimals is possible for M > M acc (R > i? acc ). In the top panel of fig. [2] 
(7 = 10~ 3 and f g — 1), planetesimal accretion proceeds in a range of i?'s in which the solid 
line (e eq ) is located below the dashed line (e acc ), that is, only if a body is larger than Ceres. 
Such large planetesimals can be formed by a different mech anism than pairwise a ccretion 



such as self-gravitational instability in turbulent eddies (e.g.. iJohansen et al.l 120071 ). When 
the disk gas is removed, accretion becomes possible for smaller planetesimals (the 2nd panel 
of fig. [2]). On the other hand, if turbulence is stronger (7 ~ 10~ 2 ), planetesimal accretion 
requires more than 1000 km-sized bodies. This appears to be an insurmountable barrier 
to accretion, even for depleted gaseous disks (f g ~ 0.1), as shown in the the 3rd panel of 
fig. [2J Finally, at large orbital radii, planetesimal accretion is even more difficult (the bottom 
panel). 

Another critical mass (radius) is the point at which ed rag (eq. [19]) and e^ in the gravity 
regime (eq. [7]) cross, 

M dis * 6 x 10-Vf (to^t) 6 - 7 (si^)" ' 6 (liu) 4 - 2 ^e, 

00/ \ -0-53 , , l/^J 

^s-SZ-do^) 2 ' 2 ^) (liu) km. 

Planetesimals with M < M,a s (R < Rdis) are disrupted by collisions down to the sizes for 
which material strength is dominant (see below). For 7 ~ 10~ 3 and f g ~ 1 (the top panel of 
fig. [2]), planetesimals with sizes larger than several km radius survive but without growing, 
while smaller planetesimals are disrupted. 

When a planetesimal is smaller than ~ 100 m in size, it is bounded by material strength 
rather than self-gravity In the regime of material strength, Qd ~ Qo(R/lcm)~ OA . The body 



0.9 
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is not disrupted if v^Qd/^k > ed rag , which is equivalent to 

M < M mat * 0.8 x 10-*V 9 (lo^)" 3 - 7 (j|-)"' 9 (Tiu) 2 ' 8M e 
i? < i? mat * 16/-"-« (^)- L25 (ai^)" ' 62 diu) 

Since in this regime, collision damping is slightly stronger than gas drag, actual values of 
M mat and R ma t are determined by \/2Qd/vk > e co ii, so they are slightly larger than the above 
estimate (see fig. 12} • When rd ra g^K ^ 1, the planetesimals' motions are coupled to that of the 
gas. The collision velocity then cannot be expressed in terms of orbital eccentricity. This 
limiting size is however much smaller than i? mat . The collision cascade would hence stop 
at M ~ M mat (R ~ -Rmat)- Regions for which the dotted lines (e^) in fig. [2] have negative 
gradients correspond to the material strength regime. In the depleted disk case, ed; s is always 
larger than e eq , so that the disruptive regions do not exist (see the 2nd panel of fig. |2J). Note 
that -R ma t oc Qq' 94 . If the planetesimals are loose aggregates so that Qq < 3 x 10 7 (the value 
for basalt rocks or water ice), the limiting size -R ma t is smaller. 

4. Accretion/destruction of planetesimals in an evolving disk 

We now put these various critical physical radii in the context of the evolution of the 
protoplanetary disk. In order to investigate the effect of departures from power-law relations 
of the surface density and temperature profiles in real disks, we also present in this section 
results obtained from a ID disk model that inc ludes an a-viscosity and photoevaporation (see 



Guillot fc Huesoll2006l ; iHueso fc Guillotl 120051 ). The parameters used in the model presented 



here are a turbulent viscosity a = 0.01 and an evaporation parameter T atm = 100 K (the 
temperature of the evaporation part of the outer disk). Another choice of the parameters 
would affect the results only marginally. 

In the numerical calculation in this section, we evaluate the equilibrium eccentricities 
e eq as a function of planetesimal radius by solving the following relation: 

r cxc = r tidal + r drag + T coll' (24) 

where the different timescales are given by eqs. (fl4"l) to (TT7|) . The survival physical radius 
for accretion i? acc is then found, for each orbital radius in the protoplanetary disk and for 
each timestep, by solving the equation, 

where e acc is given by eq. (J3J). When the mean kinetic energy is larger than the strength 
of a planetesimal, the a collision is highly erosive. We then obtain the range .R < i?dis 
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corresponding to the highly erosive collisions by solving the equation, 



/ D N A/2<5 D (i? dis ) 

e eq (-Rdis) = , (26) 

Vk 



where Qd is given by eq. (jUJ) 



Figures [3] to [5] show our results for three values of the turbulent excitation parameter, 
7 = 10~ 2 , 10~ 3 (our fiducial value), and 10~ 4 . Each figure shows, for three orbital distances, 
1, 5 and 30 AU, the planetesimal physical radii corresponding to accretive and erosive regions 
are plotted as a function of the total mass remaining in the disk. Since the disk mass 
decreases with time as a result of viscous evolution and photoevaporation, a decrease in 
disk mass corresponds to evolution in time. As shown in the previous section, planetesimal 
accretion becomes easier as the disk becomes less massive simply because the turbulent 
excitation, directly proportional to the local surface density of the gas, becomes weaker. 
However, after some point, the disk becomes too light to provide a sufficient amount of gas 
to form Jupiter-mass gas giants. 

The figures also show as thin black lines the values obtained for a disk that follows the 
slope in surface density versus orbital distance defined for the MMSN (eq. |Sj) as a function 
of a disk mass. The disk mass in this model is given by 3.4 x lO _2 / ff (r e dge/lOOAU) 1 / 2 M0, 
where r^m is the outer edge radius of the disk. Although the original MMSN model by 



Hayashil (119811 ) used r edge = 35 AU, we here adopt r e d ge = 1000 AU (for comparison, our 



fiducial alpha disk model with T atm = 100 K extends up to a maximum of 350 AU). 

These are found to be in excellent agreement with the analytical expressions derived in 
the previous section, with small differences arising from the simplifications inherent to the 
analytical approach. Larger differences are found between the power-law disk and the a-disk 
models mostly because of the difference in slopes (rflnS/rflnr = —3/2 for the former, ~ —1 
for the a-disk) which implies that a given disk mass does not correspond to the same surface 
density with two models, the difference being larger at smaller orbital distances. However, 
the qualitative features of accretive and erosive regions are similar to each other. It should 
be noted that the models also differ in their temperature profiles, but this is found to be less 
important. 

In the highly turbulent regime presented in fig. [31 a self-sustained regime of accretion 
becomes possible only when planetesimals have become very large /massive, with sizes gener- 
ally well over 100km. This case also yields a sustained area of high erosion where the average 
kinetic energies of planetesimals are above their internal energies. It is difficult to imagine 
how planetary cores can form in this context especially if they have to grow large enough to 
form giant planets. 
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With smaller perturbations from the turbulent disk (fig. HJ), planetesimals have more 
possibilities to accrete: with time, as the disk mass decreases, the inner disk rapidly moves out 
of the highly erosive regime, while erosion still remains important at large orbital distances. 
With a turbulence strength parameter 7 = 10~ 4 , corresponding to a very weak turbulence 
(fig. [5]), the presence of a highly erosive regime centered around ~ 300m planetesimals is 
limited only to the outer regions (> 30 AU), and the zone rapidly shrinks as the circumstellar 
gaseous disk disappears. 

In order to put these findings into context, we also show the mass of the disk when 
Jupiter is believed to have started accreting its gaseous envelope. These values are calculated 
by assuming that the planet growth has been limited mostly by viscous diffusion in the disk, 
with the protoplanet capturing b etween 10% and 70% of the mass flux at it s orbital distance 



in th e disk evolution model by iHueso fc Guillotl (120051 ) (for details, see iGuillot fc Hueso 



20061 ) . If giant planets have to form, at some time corresponding to the disk mass interval 
defined by the hashed areas in figs. [3] to protoplanetary cores must be already large enough 
to start accreting the surrounding hydrogen and helium gas. 

We can now define three important disk masses, and their corresponding disk ages (with 
the caution that ages are inherently model-dependent and are provided here for illustrative 
purposes only, on the basis of our particular model of a-disk evolution with photoevapora- 
tion): 

1. The maximum mass of the disk, following the collapse of the molecular cloud. This 
mass can vary quite significantly from one disk formation/evolution model to another. 
For the particular model shown here, it is of the order of 0.25 M , for an age of 0.6 Myrs. 

2. The disk mass necessary for Jupiter to grow to its present mass if it captures 10% of 
the mass flux at its orbital distance. For realistic disk models, this depends weakly 
on parameters such as a and the disk evaporation rate. In our case, it corresponds to 
M disk = 0.035 M and an age of 1.95 Myrs. 

3. The disk mass necessary for Jupiter to grow to its present mass if it captures 70% of 
the mass flux at its orbital distance. For our model, Mdi S k = 0.0054 M (about 5 times 
the mass of Jupiter) and an age of 2.85 Myrs. 

Table [T] provides the values of the physical radii that define the accretive regime and 
the highly erosive (disruptive) regime of figs. [3] to EJ namely -R acc and R^s- In our Solar 
System, the existence of Jupiter implies that either turbulence was low, the planet grew 
from a protoplanetary core formed in the inner solar system, or a mechanism was able to 
lead to the rapid formation of embryos larger than 240 km in radius at 5 AU (17 km in the 
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low-turbulence case, 1080 km in the high-turbulence case) by the time when the disk mass 
had decreased to 5 x 10~ 3 M m . In the last case, it appears that a mechanism such as the 
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standard gravitational instability (e.g., ISafronovlll969l ; iGoldreich fc Wardlll973l ) would not 



work because of the turbulence, but formation of relatively large protoplanets in eddies or 
vortexes (e.g., IJohansen et all 120071 ; iBarge fc Somerialll995l ) is a promising possibility. 



5. Conclusion and discussion 

We have investigated the critical physical radii for collisions between planetesimals to 
be accretional (R > -R acc ) or disruptive (R < -Rdis) i n turbulent disks, as functions of tur- 
bulent strength (7 ~ 0(Ap/p)), disk gas surface density, and orbital radius. The results 
presented here highlight the fact that MRI turbulence poses a great problem for the growth 
of planetesimals: generally, only those with sizes larger than a few hundred km are in a 
clearly accretive regime for a nominal value of 7 ~ 10 -3 . The others generally collide with 
velocities greater than their own surface escape velocities. For some of them, more severely 
in the kilometer-size regime, collisions are likely to be disruptive. The problem is greater 
when the disk is still massive and at large orbital distances. Also, if turbulence is stronger 
than 7 ~ 10~ 2 , planetesimal accretion becomes extremely difficult. 

However, the rate of occurr ence of extrasolar gian t planets around solar-type stars is 
inferred to be as large as ~ 20% (ICumming et al.ll2008l). and depen ds steeply on the metal- 



licity of the host star (IFischer fc Valenti 



2005 



Santos et al. 112004 ). This strongly suggests 



that the majority of extr asolar giant plan ets were formed by core accretion followed by gas 
accretion onto the cores (jlda fc Lin! 120051 ) . Thus, planetesimals should commonly grow to 
planetary masses before the disappearance of gas in protoplanetary disks. 

The possibilities to overcome the barrier are in principle as follows (their likelihood is 
commented below): 



1. Large M: Large planetesimals with sizes of 100 to 1000 km are formed directly in 
turbulent environment by a mechanism other than collisional coagulation, jumping 
over the erosive regime for physical radii. 

2. Small S g : Planetesimals start their accretion to planet-size only after the disk surface 
density of gas has declined to sufficiently small values. 

3. Small 7: Planetesimals form in MRI-inactive regions ("dead zones") of protoplanetary 
disks. 
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Concerning point 1, the first-born planetesimals with sizes larger th an R arr may be 
form ed rapidly by an efficient capture of ~meter-size boulders in vortexes (IJohansen et al. 
20071 ). Such large planetesimals may be consistent with the size distribution of asteroids 
(Morbidelli et al. 2008). Even if the first-born planetesimals are not as large, a small 
fraction of them could continue to grow larger than R acc by accreting smaller bodies, because 
accretion is not completely cut off as soon as t>di sp ^ ^esc (there is always a small possibility 
for accretion) and the large planetesimals would not be disrupted by smaller ones. This 
possibility, however, must be examined by a more detailed growth model taking into account 
the effect of fragmentation and the size distribution of planetesimals, which we neglected in 
this paper. 

Concerning point 2, we have shown that planetesimals are most fragile at early times, 
in massive disks, and at large orbital distances. We therefore suggest that the growth to- 
wards planet sizes may be delayed due to MRI turbulence, and then proceed from inside 
out: planetesimals should start accretion first close to the star, then progressively at larger 
orbital distances, as the gas surface density declines. The possibility to delay planet for- 
mation while keeping non-migrating km-size planetesimals is noteworthy because it would 
help planetary systems resisting to type-I migration: they would grow in a gas disk that 
i s less dense, and for which m igrat ion timescales may b e considerably increased. Ilda &: Lin 
(j2008al ). lAlibert et al.l (120051 ) . and iDaisaka et al.l (120061 ) showed that type-I m igration must 



be lo wered by one to two orders of magnitude from the linear calculation (ITanaka et al. 



20021 ) to provide an explanation for the existence of a population of giant planets in agree- 
ment with observations. This "late formation" scenario is consistent with the noble gas 



enrichment in Jupiter (IGuillot fc Huesd 120061 ). However, in order to form gas giants, core 
accretion and gas accretion onto the cores must proceed fast enough to capture Jupiter-mass 
amount of gas from the decaying gas disk. Once the size of the largest planetesimals exceeds 
~ 1000km, their eccentricities are damped by tidal drag and dynamical friction from small 
bodies. Most of the other small bodies may be ground into sizes smaller than 1km and 
their eccentricities could be kept very small by gas drag and collision damping. This could 
facilitate the runaway accretion of cores to become large enough (> 10M e ) for the onset 
of runaway gas accretion. This issue also has to be addressed by a detailed planetesimal 
growth model taking into account a size distrib ution. The likelihood of relatively rapid gas 
accretion without long "phase 2" is discussed by lShiraishi fc Ida! (120081 ). If planetesimal sizes 
are relatively small, gas drag damping opens up a gap in the planetesimal disk around the 
orbit of a core and truncates planetesimal accretion onto the core. The truncation of heating 
due to planetesimal bombardment enables the core to efficiently accrete disk gas. 

Concerning point 3, the MRI inactive region ("dead zones") may exist in inner disk 
regions in which the surface density is large enough to prevent cosmic and X rays from pen- 
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etrating the disk ( Gamniielll996l ; ISano et al.ll2000l ). The preservation of a dead zone can also 



contr ibute to stall type-I migration by converting it to type-II migration (JMatsumura et al. 



20071 ) or by creating a local region with a positive radial gradient of disk pressure near 



the ice line ( Ida &: Lin 2008br). Howeve r , dead zones can be eliminated by turbul ent mix- 



ing/overshoot ( Varniere fc Tagg er 2006; Turner et all 120071 ; lllgner fc Nelson 20081) . a self- 
sustaining mechanism (jlnutsuk fc Sannl 120051 ) . and dust growth (ISano et al.ll2000l ). The last 



effect comes form the fact that small dust grains are the most efficient agents for charge 
recombination. According to grain growth, the ionization of the disk and its coupling with 
the magnetic field become stronger to activate MRI turbulence. We remark that if MRI tur- 
bulence is activated, collisions are disruptive and they re-produces small grains to decrease 
the ionization degree. This self-regulation process might maintain a marginally dead state 
and keep producing small dust grains. This might be relate d with relative ch ronological 
age difference (~ 2Myr) between chondrules and CAIs (e.g., iKita et al.ll2005l ). Whether 
dead zones exist or not is one of the biggest issues in evolution of protoplanetary disks and 
planet formation. A more detailed analysis of planetesimal accretion in turbulent disks could 
impose a constraint on this issue. 

At large orbital distances (10's of AU), the existence of a highly erosive regime that lasts 
until late in the evolution of the protoplanetary disk is an important feature of this scenario. 
It shows that the entire mass of solids is highly reprocessed by collisions, in qualitative 
agreement with the paucity of presolar grains (intact remnants from the molecular cloud 
core) found in meteorites. It also prevents the growth of large planetesimals and helps to 
maintain a large population of small grains in the disks. This is in qualitative agreement with 
observations that do not indicate a significant depletion of micron- sized grains with time, 
contrary to what would be predicted in the absence of turbulence ( JDullemond fc Dominik 
20051 : iTanaka et al-lfioosh . 



In conclusion, the existence of MRI turbulence may be a threat to planetesimal accretion. 
Given the uncertainties related to these explanations, we cannot provide a definitive scenario 
for the formation of protoplanetary cores. However, it offers several promising hints to 
explain important features of planet formation as constrained by today's observations of 
protoplanetary disks, exoplanets and meteoritic samples in the Solar System. 
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Appendix 



Here we derive the r-dependence in eq. f|T3|) . If the equation of motion is scaled by 
a reference radius r% and Tki where Tki is a Keplerian period at r\, the only remaining 
non-dimensional parameter in the equations is 

64J] r 2 



(see eqs. [4], [5], [6] in lOgihara et al.1 (120071 )). Consider the equation of motion scaled by r 1 
and Tki and that scaled by T2 and Tk2- If 7 is the same and r(n) = T(r2), these two scaled- 
equations of motion are identical and evolution of eccentricity, which is a non-dimensional 
quantity, must be identical in terms of the scaled time. Note that the magnitude of excited 
e should be proportional to 7l\ Since 

r(r) = |^f-Vr(r 1 ), (28) 

S s( r i) \ r i/ 

and Ogihara et al.'s eq. (34) derived for r = 1AU is proportional to 7r(lAU), the formula for 
arbitrary r is given by replacing a year by Tk(r) and 7 by 7(r(r)/r(lAU)) in their equation. 
As a result, 

where E 9jl is S 9 at 1AU with f g = 1 (eq. [8]) and the numerical factor was corrected as 
explained in the footnote in §3.1. Assuming the simple power-law model defined by eq. (7), 

e ~ ai ^(nu)" /4 (i^) 1/2 ' (30) 
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Table 1: Limits of the accretive regime (R > -R acc ) and the highly erosive (disruptive) regime 
R < Rdis, as a function of orbital distance and disk mass. 





1 AU 


5 AU 


30 AU 


Mdisk 

[M ] 


-'•'ace -rtdis 

[km] [km] 


-f^acc -rtdis 

[km] [km] 


-f^acc 

[km] 


-Rdis 

[km] 


Fiducia 


. case: 7 = 10~ 3 ; p p = 


= 3 


gem" 3 ; Q = 3 x 10 7 ; B = ] 


L.O 





0.25 


280. 


[1.14- 


-0.061] 


1280. 


[82. - 0.00047] 


1680. 


[1000. 


- 0.00047] 


0.035 


86. 




- — 


440. 


[5.3 - 0.026] 


850. 


[200. 


- 0.0024] 


0.0054 


46. 




— 


240. 


[1.05-0.057] 


590. 


[37.- 


- 0.0069] 



Low turbulence case: 7 = 10 ; p p — 3 gem , Q = 3 x 10, B = 1.0 



0.25 


16. 


150. 


[0.30 - 0.0099] 


590. 


[22. - 0.0099] 


0.035 


3.9 


36. 





220. 


[1.07-0.053] 


0.0054 


1.7 


17. 





103. 






High turbulence case: 7 = 10 ; p p — 3 gem , Q = 3 x 10 ; B = 1.0 



0.25 

0.035 

0.0054 




[3220. 
[1460. 
[970. 



2 x 10~ 5 ] 
- 0.00012] 
0.00033] 



High material resistance: 7 = 10 ; p p — 3 gem , Q = 10 ; B = 2.0 



0.25 


280. 


[0.26 - 


-0.15] 


1280. 


[27. - 0.0013] 


1680. 


[700. 


-0.0013] 


0.035 


86. 


— 


— 


440. 


[1.6-0.068] 


850. 


[66.- 


- 0.0068] 


0.0054 


46. 


— 


— 


240. 


[0.27-0.15] 


590. 


[12. 


-0.019] 



Low material resistance^ = 10 ; p p = 3 gem ; Q = 10 ; B = 0.3 



0.25 


280. 


[11.-0.026] 


1280. 


[550. - 


- 0.00018] 


1680. 


[1440. 


- 0.00018] 


0.035 


86. 


[0.56 - 0.087] 


440. 


[41.- 


-0.011] 


850. 


[590. - 


- 0.00097] 


0.0054 


46. 





240. 


[9.2- 


- 0.024] 


590. 


[240. 


- 0.0028] 



-20- 




Fig. 1. — Evolution of the eccentricity (left panels) and of the semimajor axis (right panels) 
as a function of time, for different planetesimal masses (top to bottom). A single planetesimal 
is integrated in a turbulent disk and 5 independent runs with different random number seeds 
for the generation of turbulent density fluctuations are shown in each panel. Top panels: We 
consider no tidal and gas drag. Middle panels: Results including gas drag for planetesimals 
of M = 10~ 6 M e . Bottom panels: Results for planetesimals of M = 10~ 12 M e . In all cases, 
we assume an initial orbital distance of 1AU, f g — 1 and 7 = 0.01. 
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Fig. 2. — Equilibrium eccentricities (e eq ) as a function of the physical radius R of planetesi- 
mals (solid lines). The values of e eq is determined by the minimum of individual equilibrium 
eccentricities, e t idai, edrag, and e co ii. The critical values for accretion and destruction, e acc and 
edis, are also plotted by dashed and dotted lines, respectively. At r = 1AU, the bulk density 
p p = 3gcm -3 is assumed, while p p = lgcm~ 3 at 5 AU. The mass of the planetesimals is given 
by M = 2.1 x 10- 3 ( J R/10 3 km) 3 (p p /3gcm- 3 )M e . 
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Fig. 3.— Evolution of characteristic physical radii of planetesimals as a function of disk 
mass, at several orbital distances in the disk: 1AU (top), 5AU (middle) and 30AU (bottom). 
The solid and dot-dashed curves correspond to the boundary radius for accretion regime 
(R acc ) and to that of the highly erosive regime (-Rdis)- Two disk models have been used: 
a simple power-law model with G?lnE 5 /<ilnr = —3/2 (eq. [8]) with an outer cut-off radius 
of 1000AU (thin lines ), and an alpha-disk model with a = 0.01 and T atm = 100 K (see 
Guillot fc Huesd 120061 ) (thick lines) . The hashed region corresponds to the range of disk 
mass (equivalently, the range of time if disk evolution is given) during which Jupiter must 
start accreting hydrogen/helium gas (assuming it grabs between 10% and 70% of the disk 
mass flux at its orbital distance). In the simulations, MRI turbulence is supposed to be high, 
with 7 = 10 -2 . We also choose Q = 3 x 10 7 and B = 1 (see eq. [6]). 
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Fig. 4. — Evolution of characteristic physical radii i? acc and i?dis of planetesimals as a function 
of disk mass. The parameters and labels are the same as those in fig. [3] but for a medium 
turbulence case (7 = 1CT 3 ) (our fiducial case). 
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Fig. 5. — Evolution of characteristic physical radii i? acc and i?dis of planetesimals as a function 
of disk mass. The parameters and labels are the same as those in fig. [3] but for a weak 
turbulence case (7 = 10~ 4 ). 



